%Matlab code to generate a supplmental figure showing realizations
%of the null and alternate models in support of "Putting the
%significance of spectral peaks on the level: implications for the
%1470-yr peak in Greenland δ18O" by P. Huybers, Feb 2022.  This
%file was run on Matlab release R2021b and requires the signal
%processing toolbox. 

%Examine synthetic scenarios
addpath Toolbox;

clear;
alpha=0.01;

%%load data
addpath ../Data;
GISP2;
timeo=G(1:end,3)/1000;
Oo=G(:,2);
pl=find(timeo>=20 & timeo<=50);
O=Oo(pl);
time=timeo(pl);
mdt=mean(diff(time));
pl=find(timeo>=0 & timeo<=70);
Oo=Oo(pl);
timeo=timeo(pl);

%%interpolate to uniform intervals
ti=time(1):mdt/10:time(end);
Oi=interp1(time,O,ti); 
Oi=smoothPH(Oi,11); 
time=ti(1:10:end); 
O=Oi(1:10:end);
dt=diff(time(1:2));

%%spectral analysis 
opts.dt=dt;
opts.nw=2;
opts.mdlV=1;
dof=2*(2*opts.nw-1);
k=dof/2;
bias=psi(0,k)+log(2)-log(k*2);
conf=log(gaminv(1-alpha,k,1/k));
confb=log(gaminv(1-alpha/(length(time)/dof),k,1/k));
Vexp=psi(1,dof/2);  %expected variance of log(P)-log(Pc)
realizations=1;
ds=0.1*1/1.470;
so=1/1.470;

%%Make realizations

lab='abcdefghijklmnop';

figure(1); clf; hold on;
for rep=1:4,
    if rep==1,
        opts.p=1;
        opts.q=0;
        P(rep).opts=opts;
    end;
    if rep==2,
        opts.p=1;
        opts.q=1;
        P(rep).opts=opts;
    end;
    if rep==3,
        opts.p=3;
        opts.q=2;
        P(rep).opts=opts;
    end;
    if rep==4,
        opts.p=5;
        opts.q=5;
        P(rep).opts=opts;
    end;
    [dummy s coeffs]=specw(detrend(O),dt,opts.nw,opts);
    opts.ar=coeffs.p;
    opts.ma=coeffs.q;
    mdl = arima('MA',opts.ma,'AR',opts.ar,'Variance',opts.mdlV,'Constant',0);
    x=simulate(mdl,length(time),'NumPaths',realizations);

    xo=detrend(x);
    do=1.3*sin(2*pi*(time/1.470+rand))';
    xoa=detrend(x+do);

    [Psd s]=pmtmPH(xo,opts.dt,opts.nw,0); 
    [Psda s]=pmtmPH(xoa,opts.dt,opts.nw,0); 

    subplot(4,3,[rep*3-2 rep*3-1]); cla; hold on;
    plot(time,xoa,'r'); 
    plot(time,x,'k'); 
    xlabel('time (ka)');
    ylabel('normalized units'); 
    h4=axis; 
    text(h4(1),h4(4),['(',lab(rep*2-1),')'],'fontsize',14); 
    set(gca,'fontsize',13);

    subplot(4,3,rep*3); cla; hold on;
    h=plot(s,Psda,'r');
    h=plot(s,Psd,'k');
    set(gca,'YScale','log','XScale','log');
    axis tight;
    h4=axis;
    plot([1 1]*1/1.147,h4(3:4),'k--'); 
    xlabel('frequency (1/ky)'); 
    ylabel('power spectral density'); 
    text(h4(1),h4(4),['(',lab(rep*2),')'],'fontsize',14); 
    set(gca,'fontsize',13);
end;

return;
print -depsc fig_supp.eps;
!pstopdf fig4_supp.eps
!open fig4_supp.pdf